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We  describe  the  use  of  SIGINT,  a new  interactive  software  package  de- 
veloped by  the  National  Institute  of  Standards  and  Technology,  which 
facilitates  the  calculation  of  time  domain  frequency  stability  in  terms  of 
the  Allan  variance,  <7y(r),  or  modcTy(r)  as  a function  of  measuring  time 
from  frequency  domain  data.  Except  for  the  graphic  output,  the  code  is 
written  in  standard  FORTRAN  77  and  runs  on  AT  compatible  comput- 
ers that  have  a math  co-processor.  It  also  runs  on  many  other  systems; 
however,  calls  to  an  available  graphics  library  will  need  to  be  substituted 
for  those  that  are  included  in  this  version.  The  program  uses  either  a user 
defined  function  for  the  input  noise  or  default  functions  that  describe  the 
noise  types  commonly  found  in  oscillators,  amplifiers,  frequency  multipli- 
ers, frequency  dividers,  and  general  signal  processing  equipment.  These 
default  functions  make  it  simple  to  analyze  the  time  domain  frequency  sta- 
bility as  a function  of  measuring  bandwidth  using  realistic  first-,  second-, 
or  third-  order  low-pass  filters  or  the  simplified  infinitely  sharp  cutoff  pa- 
rameter fh.  The  default  functions  are  also  set  up  to  examine  the  effect  of 
various  servo  parameters  on  the  performance  of  a frequency  source  locked 
to  a frequency  reference. 

Key  words:  Allan  variance;  frequency  lock  loop  analysis;  modified  Allan 
vaxictnce;  pha^e  lock  loop  analysis;  spectral  density  of  frequency  stability; 
time  domain  frequency  stability. 


1 Introduction 


SIGINT  is  an  interactive  software  package  designed  to  facilitate  the  calculation  of 
time  domain  frequency  stability  from  data  in  the  frequency  domain.  Both  the  AIIeui 
variance  (two-sample  variance),  <Ty(7"),  and  the  modified  Allan  Vciriance,  moday(r), 
can  be  calculated  as  functions  of  measurement  time.  Except  for  the  graphic  output, 
the  code  is  written  in  standard  FORTRAN  77  to  run  on  AT  and  compatibles  that 
have  a math  co-processor  (see  §4  for  details  on  portability).  The  actual  output  of  the 
software  is  the  square  root  of  the  variance.  The  equations  are  [l,2,3,4]' 


o-v(nro)  = 


SM) 


sm*(irfnTo) 

{nfriToy 


2 

2 


,and 


modcTy  (nro) 


~rr~2  (” / sm*{vfnTo)  df 

U^TT^Tq  [ •'o  P 

ffh^  sjf) 

+ 2 y X]  (^  - k) — 7^  cos(27r//:ro)  sin^(7r/nro)  df 
° jk=l  J 


where  the  measurement  time  is  given  by  tito  with  ro  being  the  minimum  measurement 
time. 


The  input  parameters  are  expressed  in  terms  of  the  spectral  density  of  frequency 
fiuctuations,  Sy{f).  Conversion  from  phase  noise  to  Sy{f)  is  simply 

Sy{f)  = ^S,{f),and 

sAf)  = jMf)- 

You  may  choose  to  supply  your  own  function  Sy(f),  or  to  use  the  functions  which 
are  built  into  the  package.  Use  the  parameter  SELSY  (see  below)  to  communicate 
this  choice  to  the  package.  If  you  use  your  own  function,  you  must  modify  the  empty 
FORTRAN  function,  SYF(f),  within  the  package.  The  built  in  function  has  the  form: 

Cif~^  + ^2/“^  + C3  + C4/  + c^p 


Sy{f)  = 


1 


The  values  of  the  C,-  are  stored  in  array  C.  The  filter  function  K{f)  must  be  chosen 
from  the  following  four  functions  by  setting  the  value  of  SELK  (SELK=0,  1,  2 and  3 
respectively) . 


K(f)  = 1, 

K(f)  = + > 


The  coefficients  Ki  must  be  input  in  the  array  CK. 

M[f)  = (1+Mi/(1  + M2/)(1+M3/))' 

where  any  or  all  of  the  M,  can  be  0.  These  M,  are  stored  in  the  CM  array. 

This  accommodates  the  most  conunonly  encountered  types  of  random  noise  found  in 
oscillators,  amplifiers,  frequency  multipliers,  frequency  synthesizers  and  general  signal 
processing  equipment.  In  addition  the  upper  cutoff  frequency  of  the  integral  can  be 
treated  as  infinitely  sharp  by  setting  fh  equal  to  the  equivalent  noise  bandwidth  of 
the  time  domain  configuration  simulated  by  the  calculations,  or  it  can  be  treated 
more  realistically  using  the  appropriate  first-,  second-,  or  third-order  low-pass  filter 
function.  This  option  is  implemented  by  choosing  the  appropriate  constants  for  Mi, 
M2,  and  M3,  and  setting  fh,  to  a sufficiently  large  value  that  the  termination  of  the 
integration  at  this  value  does  not  significantly  bias  the  results. 

Servo  analysis,  such  as  the  examination  of  the  effect  of  various  servo  parameters  on 
the  locking  of  an  oscillator  to  a frequency  reference,  is  facilitated  by  choosing  the 
appropriate  form  of  K(f).  The  noise  type  assumed  for  the  reference  source  is  white 
frequency  modulation.  An  example  of  this  is  given  in  the  Appendix.  In  some  special 
servo  cases,  it  may  be  necessary  to  divide  Sy[f)  into  two  segments  and  run  them  as 
separate  cases  , or  you  can  enter  your  own  form  for  Sy(f). 

The  rest  of  the  pcirameters  are  described  in  the  Input  Parameter  section,  below. 
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2 Running  the  Program 


The  SIGINT  package  of  routines  is  set  up  for  interactive  use.  It  is  machine  independ- 
ent with  the  exceptions  noted  in  the  Portability  section  below.  The  package  will  first 
display  the  values  of  the  parameters  that  define  the  default  case.  It  then  presents  a 
menu  which  asks  whether  you  want  to  run  that  case,  be  prompted  to  input  an  entire 
new  set  of  parameters,  or  modify  directly  the  current  values  of  some  parameters.  This 
will  be  referred  to  later  as  the  “input  menu” . 


2.1  Input  Modes 

In  either  of  the  modes  described  below,  values  may  be  entered  in  either  fixed-  or 
fioating-point  format.  If  fioating-point  is  used,  the  decimal  point  is  necessary;  for 
example,  you  can  enter  2.e-3  but  not  2e-3.  In  some  cases,  when  an  integer  is  being 
entered,  the  decimal  point  must  be  omitted.  In  all  cases,  the  program  recovers  from 
incorrect  input,  outputs  a message  attempting  to  diagnose  the  problem,  and  allows 
re-entry  of  the  value. 


2.1.1  The  “prompting”  input  mode 

If  you  elect  to  be  prompted  for  the  values  of  the  variables,  the  program  will  lead  you 
through  a series  of  questions  and  answers  by  which  a new  case  will  be  defined.  An 
advantage  of  this  mode  is  that  it  is  not  necessary  to  know  the  names  of  the  variables 
which  the  program  uses. 


2.1.2  The  “direct”  input  mode 

K,  instead,  the  “direct”  mode  of  input  is  chosen,  values  axe  changed  by  typing  one 
or  more  lines  containing  variable  names  followed  by  “=”  and  the  new  values  of  the 
variable.  Different  entries  can  be  separated  by  a comma  or  blanks.  The  last  line  must 
be  terminated  by  a “$”  or  a For  example,  you  might  type  the  following  lines: 


NRANGE=2,  NLOW  = 10 
SELSY=1  CK=1.,2.,3.  $ 


This  mode  will  usually  be  faster  than  the  “prompting”  mode. 

Notice  that  if  you  want  to  get  out  of  this  mode  without  changing  anything,  you  can 
simply  type  a “$”  or  a . 

When  you  finish  with  either  mode  of  input,  the  current  set  of  parameters  will  be 
displayed  and  you  will  be  returned  to  the  input  menu.  At  this  point  you  can  either 
run  the  Ccise  or  further  modify  the  parameters  by  again  entering  one  of  the  input 
modes.  You  can  continue  to  revise  and  examine  the  parameters  until  you  are  satisfied 
with  them  and  then  run  the  case. 


2.2  Output  Modes 

When  the  integration  is  completed,  you  will  be  presented  with  another  menu  which 
will  be  called  the  “output  menu”.  This  menu  asks  whether  to  print  results,  plot  them, 
compute  another  integral,  or  quit. 

As  the  program  computes  an  integral,  it  saves  the  defining  parameters  and  the  results 
in  internal  data  structures.  The  results  are  in  the  form  of  a table  of  values  of  r 
(r  = titq)  and  corresponding  values  of  the  square  root  of  the  integral  (variance).  So 
eaeh  time  you  come  to  the  output  menu,  the  results  of  all  of  the  cases  you  have  done 
during  the  run  are  available  to  be  printed  and/or  plotted. 


2.2.1  Printing 

On  AT  compatible  computers  you  can  either  print  to  the  screen  or  to  your  printer 
(see  Portability  section  for  behavior  on  other  systems) . 


2.2.2  Plotting 

Plotting  will  work  only  on  AT  compatible  computers.  On  each  frame  it  is  possible  to 
plot  either  one  or  several  cases.  A hard  copy  of  the  plot  Ccin  be  obtained  if  you  have 
a memory-resident  utility  that  allows  you  to  use  the  print  screen  key  to  transfer  the 
graphics  screen  image  to  your  printer.  The  plot  package  will  plot  directly  to  certain 
printers.  Examples  of  the  plots  are  found  in  Figures  1-7.  Additional  details  can  be 
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obtained  by  contacting  the  Applied  and  Computational  Mathematics  Division,  NIST, 
325  Broadway,  Boulder,  Colorado  80303-3328.  Attention:  Abbie  O’Gallagher. 


2.2.3  Quitting  and  Saving  Results 

When  you  choose  to  quit,  you  will  be  asked  if  you  want  to  save  the  results.  If  you 
say  yes,  the  defining  parameters,  along  with  the  results  of  each  of  the  cases  you  have 
done  during  that  nm,  will  be  saved  in  a disk  file  called  SAVOUT  on  the  default  disk 
and  directory. 


2.3  Input  Parameters 

The  following  is  a list  of  the  input  variables  that  must  be  set  in  order  to  compute 
an  integral.  The  integral  can  be  computed  for  a single  value  of  r (r  = nro),  or  for  a 
sequence  of  values  of  r.  In  the  latter  case,  the  integral  can  be  plotted  as  a function 
of  T. 


INTGRAL: 


Integer  which  selects  the  integral  to  be  computed.  INTGRAL  = 1 for 
the  a integral;  INTGRAL  = 2 for  the  mod<7  integral.  The  default  value 


is  INTGRAL  = 1. 


NRANGE:  Integer  which  determines  the  range  of  values  of  the  parameter  n for 

which  the  integral  is  computed.  Set  NRANGE  = 1 to  compute  the 
integral  for  a single  value  of  n (n  = NLOW).  Set  NRANGE  ==  2 to 
compute  the  integral  for  the  sequence  of  values  obtained  by  starting 
with  NLOW  and  doubling  the  current  value  to  obtain  the  next  one  until 
reaching  the  last  such  value  not  exceeding  NHIGH.  Thus  if  NRANGE 
= 2,  NLOW  = 1,  and  NHIGH  = 20,  then  the  integral  will  be  computed 
for  n=(l,  2,  4,  8,  16).  For  NRANGE  = 3,  the  integral  will  be  evaluated 
by  decades  with  five  points  per  decade.  Thus  if  NRANGE  = 3,  NLOW 
= 1,  NHIGH  100,  then  the  evaluation  will  be  done  for  n = (1,  2,  3, 
5,  7,  10,  20,  30,  50,  70,  100).  The  default  value  is  NRANGE  = 3. 

NLOW:  The  lower  limit  of  the  range  of  values  of  n for  which  the  integral  is 

computed  (see  the  discussion  of  NRANGE).  The  default  value  is  NLOW 

= 1. 
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NHIGH:  The  upper  limit  of  the  range  of  values  of  n for  which  the  integral  is 

computed.  The  default  value  is  NHIGH  = 1000. 

FH:  The  upper  limit  of  the  integral,  fh.  The  default  value  is  FH  = 3. 

SELSY:  Integer  to  select  the  function  Sy[f).  If  SELSY  = 1,  you  must  supply 

a double-precision  function  named  SYF(f),  with  double-precision  argu- 
ment f.  (You  must  replace  the  empty  FORTRAN  function  SYF(f)  that 
is  provided.)  If  SELSY  = 2,  then  you  must  select  one  of  four  built-in 
functions  by  giving  values  of  the  input  parameter  SELK,  and  the  arrays 
C,  CK,  and  CM  as  given  in  the  discussion  of  Sy(f)  above.  The  default 
value  of  SELSY  is  SELSY  = 2. 

SELK:  Integer  parameter  which  selects  the  function  K(f).  See  discussion  above. 

The  default  value  is  SELK  = 3. 

C:  A double-precision  array  of  dimension  C(6)  containing  the  coefficients 

of  the  function  Sy{f).  The  default  values  axe  C = 2.e-24,  0,  0,  0,  0,  0. 

CK:  A double-precision  array  of  dimension  CK(3)  containing  the  coefficients 

of  the  function  K(f).  The  default  values  axe  CK  = 10,  40,  100. 

CM:  A double-precision  array  of  dimension  CM(3)  containing  the  coefficients 

of  the  function  M(f).  The  default  values  axe  CM  = 0 ,0,  0. 

TAUO:  A double-precision  scalar  parameter  of  the  integral.  The  default  value 

is  TAUO  = 1. 


2.4  Internal  Constants 


It  is  not  necessciry  to  be  aware  of  the  constants  described  in  this  section  unless  you 
want  to  alter  the  behavior  of  the  code.  In  that  case  it  may  be  necessary  to  change  the 
values  to  which  these  constants  are  set.  This  may  be  done  by  editing  the  FORTRAN 
to  change  the  assignment  statements  that  set  these  values.  Recompiling  will  also  be 
necessary.  The  pertinent  assignment  statements  can  be  found  at  the  beginning  of  the 
main  program. 
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EPSREL:  The  relative  error  tolerance  for  QUADPACK  routines.  It  is  set  at 

2 * 10"^. 

EPSTRN:  The  truncation  error  tolerance  used  by  the  MDSIGY  routine  to  deter- 

mine, for  a given  value  of  the  variable  of  integration,  whether  the  esti- 
mate of  the  integral  over  the  remainder  of  the  interval  is  small  enough 
to  indicate  that  the  remainder  can  be  ignored.  See  discussion  of  the 
numerical  method,  below.  EPSTRN  is  set  at  2 * 10“^. 


NCYCLE:  The  number  of  cycles  of  the  numerator  of  Sy{f)  on  each  side  of  a sin- 

gularity which  are  evaluated  by  the  non-oscillatory  quadrature  routine 
(see  section  3).  NCYCLE  is  set  to  8. 


IDB:  Switch  to  control  debug  printout  from  the  routines  SIGMAY  and  MD- 

SIGY. It  is  set  to  0.  Set  it  to  1 to  get  debug  printing. 


3 The  Numerical  Method 


The  sigma  integral. 


(^y{nro)  = 


2 


sin^(7r/nro) 

(7r/nro)2 


I 

3 


1 


can  be  treated  as  an  oscillatory  integral  except  in  the  neighborhood  of  the  origin. 
Therefore,  the  integral  is  computed  by  using  one  approximation  in  the  neighborhood 
of  the  origin  and  a second  method  away  from  the  origin;  that  is, 

where 


h 

h 

Piif) 

m) 


•'o 

/ ^2  (/)  sin^  (vfriTo)  df, 

a 


Sy{f) 


Sy{f) 


sin^(7T  friTo) 
(TrfriToy 


(7r/nro)2’ 


,and 
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with  a chosen  so  that  the  first  integral  is  tciken  over  NCYCLE/2  cycles  of  the  sine 
function,  provided  r = uto  is  large  enough.  Otherwise  it  is  taken  over  2/3  of  the 
interval  fh.  That  is,  a = mimmum(2/fc/3, NCYCLE/r).  E the  function  Fi  grows 
more  slowly  than  f~^  at  / = 0,  then  the  first  integral  exists.  We  assume  that  Sy(f) 
grows  no  more  rapidly  than  /“^,  and  thus  the  integrand  in  this  first  integral  has  a 
removable  singularity  and  can  be  computed  using  a general  purpose  quadrature  rou- 
tine from  the  subroutine  library.  We  use  the  DQAGE  routine  from  the  QUADPACK 
[7]  library. 

The  second  integral  is  oscillatory.  It  can  be  transformed  into  a canonical  form  by  use 
of  the  trigonometric  substitution, 

3 1 1 

sin'^(7r/r)  = cos(2tt fr)  + - cos(47r/r), 

8 2 8 

and  can  therefore  be  written  as  the  sum  of  three  integrals, 

h = hi  + h2  + hz‘ 

Since  a is  chosen  large  enough  to  avoid  the  singularity  at  the  origin,  the  first  of 
these  three  integrals  can  be  computed  with  the  standard  quadrature  routine  DQAGE. 
The  remaining  two  integrals  are  oscillatory  and  are  therefore  computed  using  the 
DQAWOE  routine,  also  from  QUADPACK.  DQAWOE  is  designed  for  integrals  of 
the  form  / f{x)  sin(o;x)dr  or  / f(x)  cos(ux)dx  ,where  u can  be  very  large.  These 
sigma  integrals  may  have  values  of  nro  which  easily  exceed  10^  (here  uto  corresponds 
to  cj). 

In  all  cases  we  ask  that  the  QUADPACK  routines  compute  the  integrals  with  a 
relative  error  less  than  2 * 10“^.  This  may  be  wasteful,  since  one  of  these  integrals 
may  dominate  the  remainder,  however  we  could  not  be  certain  of  the  dominance  in 
all  cases. 

Extensive  tests  of  the  numerical  calculation  for  CTy(r)  have  been  performed  for  all  the 
common  noise  types,  Ci  — C5  [1,2,3].  The  results  agree  with  the  asymptotic  form 
(27r/^r  ':$>  1)  to  better  than  11%  as  long  as  2nfhT  is  larger  than  6.3  and  0.1%  for 
27t fhT  > 62.8.  We  believe  that  the  numerical  technique  yields  more  accurate  results 
since  the  analytical  calcuations  are  valid  only  for  2'KfhT  ^ 1.  The  asymptotic  forms 
are  considered  in  more  detail  in  Figures  1-7  and  Tables  I and  II. 
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The  mod  sigma  integral, 


moday(nro) 


= -rV?  ( df 

U^TT^Tq  ( •'O  P 

+ 2 (^  “ cos{2nfkTo)  sm^{7rfnro)  df 

° jb=i  J 
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is  much  more  difficult  to  approximate  because  the  sum  over  k includes  many  terms 
when  n is  large.  The  sum  over  k is  eliminated  by  using  the  following  identity, 


MO,  ” Ism^^ 

Z(n-fc)cos0fc  = -• ; + 7-5-2T- 
k=i  2 2 sm  j 


This  identity  is  derived  by  summing  a geometric  series  composed  of  complex  expone- 
tials.  Since  the  real  part  of  a complex  exponential  is  the  cosine  term,  this  series 
reduces  to  the  expression  on  the  left,  above.  This  can  be  summed  explicitly,  yielding 
the  final  formula.  Application  of  this  identity  to  the  original  integral  leads  to  the 
following  form  of  the  integral: 


modcTy(nro) 


2 ff>‘  Sy(f)  sin^(7rron/)  1 “ 
n*n^T^  ■^0  P sm^{iTTof) 


With  this  simplification  mod<jy(nro)  can  be  broken  up  and  evaluated  in  a way  that 
is  analogous  to  that  described  above  for  the  <7y(nro)  integral.  However  mod(7y(nro) 
presents  the  additional  problem  that  it  has  a singularity  for  each  integral  value  of 
Tof  instead  of  just  one  singularity  at  / = 0.  It  is  therefore  necessary  to  break  up  the 
interval  of  integration  into  p * tq  subintervals  and  sum  the  values  of  the  integral  over 
each  of  these  smaller  intervals.  If  is  large,  this  is  a large  number  of  subintervals.  But 
in  many  cases,  the  later  subintervals  contribute  negligibly  to  the  value  of  the  integral. 
For  this  reason  an  estimate  of  the  error  produced  by  neglecting  these  contributions 
is  made  every  few  subintervals.  If  the  ratio  of  this  error  to  the  computed  value  of  the 
integral  up  to  that  point  is  small  (less  than  EPSTRN),  then  no  further  calculation  is 
done  and  the  value  returned  by  the  program  as  modcTy  (nro)  is  the  square  root  of  the 
sum  of  the  integrations  done  over  the  subintervals  up  to  that  point. 

Extensive  tests  of  the  numerical  calcuation  for  mod(7y(r)  have  been  carried  out  for 
the  common  types  of  noise  (Ci  — Cs)  and  compared  to  the  results  of  [3]  , obtained 
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analytically  and  the  results  of  [4]  obtained  by  computer  simulated  noise.  This  has  been 
done  by  comparing  R(n)  the  ratio  of  moda^ (r)  to  as  a function  of  the  number  of 
samples  averaged  together  for  mod<7y(r),  obtained  by  the  three  methods.  The  results 
obtained  using  SIGINT  for  R(n)  are  shown  in  figure  1 and  table  1.  This  agrees  exactly 
with  the  results  of  both  [3]  and  [4]  for  a = —2,  a = 0 and  a = 2 (ignoring  the  obvious 
typographical  errors).  The  results  for  a = 1 agree  for  27rfhTo  = 10'*  - the  only  case 
addressed  by  [3].  Also  shown  are  our  results  for  a = 1 and  27rfhTo  = 3,  10,  100.  Our 
results  differ  considerably  from  [3]  for  a = — 1,  but  agree  with  the  asymptotic  value 
of  [4].  Since  we  have  obtained  excellent  agreement  for  so  many  other  cases  using  the 
same  numerical  code,  we  believe  that  the  results  presented  in  figure  1 cuid  table  1 are 
more  accurate  than  the  earlier  work. 


4 Portability 

The  SIGINT  code  is  standard  FORTRAN  77  and  is  portable  with  the  following  ex- 
ceptions: 

• The  graphics  library  used  is  David  Kahaner’s*  GRAPH.LIB,  which  runs  only 
on  PC  compatible  computers  with  one  of  the  following  display  adapters: 

— IBM  Color  Graphics  Adapter^ 

— IBM  Enhanced  Graphics  Adapter 

— IBM  Professional  Graphics  Adapter  in  CGA  emulation  mode 
— IBM  PS/2  Video  Gate  Array  in  EGA  emulation  mode 
— Hercules 
- Video-7 

When  moving  the  code  to  another  system,  it  will  be  necessary  either  to  “com- 
ment out”  the  calls  to  the  plot  routines  and  run  the  code  without  graphic  output 

^Scientific  Computing  Division,  National  Institute  of  Standards  and  Technology,  Gaithersburg, 
MD  20899 

^Identifacation  of  some  commercial  materials  has  been  necessary  in  this  report.  In  no  case  does 
such  identification  imply  recommendation  or  endorsement  by  the  National  Institute  of  Standards  and 
Technology,  nor  does  it  imply  that  the  materials  are  necessarily  the  best  available  for  the  purpose. 
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or  to  change  those  graphics  calls,  all  of  which  are  in  SUBROUTINE  PLOT,  to 
call  an  available  graphics  library. 

• On  any  system,  numerical  results  can  be  printed  on  the  screen  but  if  the  option 
of  using  the  printer  is  chosen  on  a system  other  than  an  AT  compatible,  the 
results  will  instead  be  written  to  a file  called  “PRN”.  (“PRN”  is  the  name  DOS 
uses  for  the  printer.)  Results  will  not  accumulate  in  this  file.  Only  the  last 
printing  will  be  there  when  the  program  terminates. 

• The  code  is  written  in  double  precision,  which  may  be  unnecessary  on  machines 
with  larger  word  sizes. 

• On  UNIX  systems,  carriage  control  characters  may  appear  on  the  screen  instead 
of  affecting  the  output  as  intended. 

• The  VAX  reports  underfiow  conditions  but  terminates  normally  and  the  results 
are  imaffected. 


The  entire  code  has  been  run  extensively  on  AT  compatible  computers  using  RM/FOR- 
TRAN.  In  addition,  except  for  the  graphic  output,  it  has  been  tested  on  a VAX 
11/785,  SUN  3/180,  MASSCOMP  MC5500-PEP  and  Cyber  840  and  855.  It  is  ex- 
pected that  it  will  run  on  most  other  systems.  If  problems  are  encountered  in  moving 
it  to  another  system,  an  attempt  will  be  made,  subject  to  the  availability  of  resources, 
to  overcome  them.  For  such  assistance,  please  contact  the  Applied  and  Computa- 
tional Mathematics  Division,  NIST,  325  Broadway,  Boulder,  Colorado  80303-3328. 
Attention:  Abbie  O’Gallagher. 
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Figure  1.  cry(r)  versus  r for  the  five  common  power-law  noise  types  in 
the  limit  that  27rf^  is  large  compared  to  1 and  an  infinitely  sharp 
filter  is  used.  Curve  a is  for  random-walk  frequency  modulation,  Sy(f) 
= h_2  f~^.  Curve  b is  for  flicker  frequency  modulation,  Sy(f)  = h_i  f" 
Curve  c is  for  white  frequency  modulation,  Sy(f)  = hQ.  Curve  d is 
for  flicker  phase  modulation,  Sy(f)  = h]^  f.  Curve  e is  for  white  phase 
modulation,  Sy(f)  = h2  f^. 
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Figure  2.  ay(r)  for  white  phase  modulation  (a  = 2) 
measurement  time,  r,  and  measurement  bandwidth, 
have  an  infinitely  sharp  filter  with  width,  fj^  = 16 
fh  = 0.0016  Hz  respectively.  Curves  d and  e have  a 
width,  fj^  = 0.016  Hz  and  f^  = 0.0016  respectively. 


i I I I I I I 

l.E+03 


as  a function  of 
Curves  a,  b,  and  c 
Hz,  fh  = 0.016  Hz, 
single  pole  filter 
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Figure  3.  Oy{r)  for  flicker  phase  frequency  modulation  (a  = 1)  as  a 
function  of  measurement  time,  r,  and  measurement  bandwidth,  f]^.  Curves 
a,  b,  and  c have  an  infinitely  sharp  filter  with  width,  fj^  = 16  Hz,  f^ 
= 0.016  Hz,  fj^  = 0.0016  Hz  respectively.  Curves  d and  e have  a single 
pole  filter  width,  f^  = 0.016  Hz  and 
f^  = 0.0016  respectively. 
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Figure  4.  Oy{r)  for  white  frequency  modulation  (a  =0)  as  a function 
of  measurement  time,  r,  and  measurement  bandwidth,  f^.  Curves  a,  b, 
and  c have  an  infinitely  sharp  filter  with  width,  f^  = 16  Hz,  fh  = 
0.016  Hz,  fh  = 0.0016  Hz  respectively.  Curves  d and  e have  a single 
pole  filter  width,  fh  = 0.016  Hz  and  fh  = 0.0016  respectively. 
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Figure  5.  Oy(r)  for  flicker  frequency  modulation  (a  = -1)  as  a 
function  of  measurement  time,  r,  and  measurement  bandwidth,  f]-^.  Curves 
a,  b,  and  c have  an  infinitely  sharp  filter  with  width,  f^  = 16  Hz,  f^ 
= 0.016  Hz,  fh  = 0.0016  Hz  respectively.  Curves  d and  e have  a single 
pole  filter  width,  fj^  = 0.016  Hz  and  f^  = 0.0016  respectively. 
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Figure  6.  Oy{r)  for  random  walk  frequency  modulation  (a  = -2)  as  a 
function  of  measurement  time  and  measurement  bandwidth,  fj-j,  for  tt. 
Curves  a,  b,  and  c have  an  infinitely  sharp  filter  with  width,  f]-j  = 16 
Hz,  f^  = 0.016  Hz,  f^  = 0.0016  Hz  respectively.  Curves  d and  e have  a 
single  pole  filter  width,  f^  = 0.016  Hz  and  fj^  = 0.0016  respectively. 
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R(n) 


n,  Number  of  Samples  Averaged 


Figure  7.  Ratio  of  mod(7y(r)^  to  ay(r)^  as  a function  of  n,  the  number 
of  points  averaged  to  obtain  Moday(r) . The  measurement  time  r = nrQ, 
where  tq  is  the  minimum  data  interval. 
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TABLE  I.  Ratio  of  modCTy(r)  to  0^(7)  vs  n,  for  common  power- law  noise  types 
Sy(f)  = h^f®.  n is  the  number  of  time  or  phase  samples  averaged  to  obtain 
moday(r  = nrg ) where  Tq  is  the  minimum  sample  time,  and  is  2n  times  the 
measurement  bandwidth  fj^  . 


modo„(x) 

R(n)=  ^ ^ = ^^0 

o;(t) 


n 

a = -2 

a = - 1 

a = 0 

= 3 

a 

= 10 

= +1 

= 100 

0 

It 

0 

a = +2 

1 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

2 

0.859 

0.738 

0.616 

0.568 

0.543 

0.525 

0.504 

0.500 

3 

0.840 

0.701 

0.551 

0.481 

0.418 

0.384 

0.355 

0.330 

4 

0.831 

0.681 

0.530 

0.405 

0.359 

0.317 

0.284 

0.250 

5 

0.830 

0.684 

0.517 

0.386 

0.324 

0.279 

0.241 

0.200 

6 

0.828 

0.681 

0.514 

0.349 

0.301 

0.251 

0.214 

0.167 

7 

0.827 

0.679 

0.507 

0.343 

0.283 

0.235 

0.195 

0.143 

8 

0.827 

0.678 

0.506 

0.319 

0.271 

0.219 

0.180 

0.125 

10 

0.826 

0.677 

0.504 

0.299 

0.253 

0.203 

0.160 

0.100 

14 

0.826 

0.675 

0.502 

0.274 

0.230 

0.179 

0.137 

0.0714 

20 

0.825 

0.675 

0.501 

0.253 

0.210 

0.163 

0.119 

0.0500 

30 

0.825 

0.675 

0.500 

0.233 

0.194 

0.148 

0.106 

0.0333 

50 

0.825 

0.675 

0.500 

0.210 

0.176 

0.134 

0.0938 

0.0200 

100  0.825 

0.675 

0.500 

0.186 

0.159 

0.121 

0.0837 

0.0100 

LiBic  0.825 

0.675 

0.500 

r 

3.37  1 

1/n 

1.04  +3  T j 
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TABLE  II.  Asymptotic  forms  of  crj(r)  for  various  power-iaw  noise  types  and  two 
filter  types.  Notei  (^/27t  = fj^  is  tlie  measurement  system  bandwidth  often 
called  the  high-frequency  cutoff.  In  = loge. 


Nane  of  Noise 

a 

s,(f) 

o^ir) 

Whr»l 

Infinite  Sharp 
Filter 

Whr»l 

Single  Pole 
Filter 

i^r«i 

Infinite  Sharp 
Filter 

WhT«l 

Single  Pole 
Filter 

White  Phase 

2 

3fhh2 

K 

3fhh2 

2/5^2  fi,5r2h2 

fh"h2 

(2jr)2r2 

(2>r)2r2 

2r 

Flicker  Phase 

1 

(1.038  + 31n(whT))hi 

(31n(whT 

))hi  hw^fi^^r^hi 

2fi,2(ln(2))hi 

hi  t 

(2»r)2r2 

(2rr)2 

r2 

White  Frequency 

0 

ho 

ho  

2r 

ho 

2r 

2/3jr2fj^2^h(j 

Flicker  Frequency 

-1 

h.if-'  2(ln(2))h.i 

2(ln(2))h 

^2f,2r2h.i 

8jt^  r^h.  1 

Random-Walk 

Frequency 

-2 

2n^ rh. 2 

H -F"  ^ 

rh-2 

2n^  fj,  T^h.  2 

2jt^  fj,  r^h.  2 

2 ^ 

3 

3 

i 
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Appendix 


Stability  of  Frequency  Locked  Loops 
Fred  L.  Walls 

National  Institute  of  Standards  and  Technology 
Boulder,  Colorado  80303 

1 ■ Introduction 


Passive  frequency  standards  are  characterized  by  the  use  of  a reference 
resonance  to  stabilize  the  frequency  of  an  external  probe  oscillator.  A 
common  configuration  is  shown  in  Fig.  1 [1-6].  The  probe  oscillator 

generally  has  phase  modulation  imposed  on  the  carrier  in  order  to 
interrogate  the  resonance  with  a minimum  of  offset.  The  resulting 
amplitude  modulation  is  demodulated  to  yield  an  error  curve  that  is 
essentially  the  derivative  of  the  resonance.  Although  Fig.  1 shows  a 
transmission  system,  similar  schemes  are  sometimes  possible  in  reflection 
[5].  The  error  signal  from  the  demodulator  is  used  to  steer  the  probe 
signal  toward  the  center  of  the  resonance  line  [1-9].  For  analysis  times 
longer  than  one  period  of  the  modulation  cycle,  and  under  the  condition 
that  the  probe  oscillator  wanders  less  than  the  half  width  of  the  error 
curve  in  the  loop  attack  time,  we  can  treat  this  curve  as  approximately 
static[l-4].  Near  line  center  the  loop  error  voltage  at  the 

synchronous  detector,  is  approximately  = k (1^0°  - where  k is 

the  slope  of  the  error  curve,  is  the  resonance  frequency  of  the 

reference,  and  is  the  open  loop  frequency  of  the  probe  oscillator  and 

is  the  detector  noise.  If  we  now  close  the  loop  with  gain  G(f),  it 
can  be  shown  that  the  spectral  density  of  fractional  frequency 
fluctuations  Sy(f)  for  the  probe  source  becomes 


MODULATION 

REF 


Figure  1 . Generalized  block  diagram  of  a probe  source  locked  to  a 
reference  resonance. 


Contribution  of  the  U.S.  Government;  not  subject  to  copyright. 
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+ s/(f) 


(1) 


Sy(f) 


G(f)  ^ 

2 

L R 

C 1 

lG(f)+lJ 

G2(f)  J 

where  Sy° (f)  is  the  open-loop  spectral  density  of  fractional  frequency 
fluctuations  of  the  probe  source,  Sy^(f)  is  that  of  the  reference,  and 
Sy^(f)  is  that  of  the  detector  and  interrogation  noise  referred  to  the 
demodulator  output.  Cutler  [10]  has  pointed  out  that  noise  in  the  local 
oscillator  at  the  2nd  harmonic  of  the  modulation  frequency,  which  is 
usually  ignored,  causes  a time  varying  frequency  offset  that  is 
undistinguishabie  from  reference  noise.  This  sets  the  lower  limit  to  the 
interrogation  noise  and  often  sets  the  lower  limit  of  the  noise 
performance  of  the  local  oscillator  necessary  not  to  degrade  the  overall 
performance.  The  magnitude  of  the  2nd  harmonic  noise  modulation  in 

radians/s  is  estimated  in  appendix  B of  [1]  to  be  kg  = 27rr/(  Sy  (^/tt)  ) ^ ^ ^ , 
where  n/(27r)  is  the  modulation  frequency.  This  leads  to  an  interrogation 
noise  term  which  is  of  order  [1,4,10]  Sy^(f)  = 1/(167t)  Sy°  (fi/Tr)  . For 
large  values  of  G(f)  , Sy(f)  of  the  probe  source  reflects  that  of  the 
reference  plus  the  added  noise  of  the  detection  system  [1-6]. 

The  primary  goal  of  this  paper  is  to  investigate  the  effect  of  various 
realistic  forms  of  G(f)  on  the  spectral  density  of  frequency  and 
fractional- frequency  stability.  It  will  be  shown  that  mod  <Jy(r)  [11,12] 
is  better  suited  than  the  traditional  two-sample  or  Allan  Variance 
[8],  for  evaluating  the  locked  performance  when  the  frequency  stability 
of  the  reference  is  much  greater  than  that  of  the  local  oscillator 
[7,8,11,12]. 

2 . The  Effect  of  Different  Forms  of  Servo  Gain 


Figure  2 shows  the  effect  of  locking  a probe  source  with  Sy°(f)  = 2 x 
10"^®/f^  + 1 X 10“  ^^/f  + 2 X 10“  ^°f^  (which  roughly  corresponds  to  that 
of  a low-noise  5 MHz  quartz  oscillator)  to  a reference  resonance  with 
Sy^(f)  of  2 X 10“^°  and  gy^(f)  = 5.6  x 10“^®.  Sy^(f)  is  the  estimated 
interrogation  noise  for  a modulation  frequency  of  47  Hz  [1,4,10].  Curves 
A,  B,  and  C show  the  effect  of  using  a first-,  second-,  or  third-  order 
loop  each  having  a loop  bandwidth  of  approximately  0.1  Hz  or  an  attack 
time  of  1.6  s[l].  The  solid  Curves  A,  B,  and  C of  Fig.  3 show  o^{t) 
calculated  for  curves  A,  B,  and  C from  Fig.  2 with  a noise  bandwidth  of  3 
Hz,  The  improvement  in  stability  of  the  local  oscillator  due  to  the 
servo  scales  roughly  as  r/r^  where  is  the  attack  time  and  r is  the 
measurement  time.  Considerable  insight  into  the  effect  of  various  gain 
stages  on  the  frequency  stability  can  be  obtained  by  considering 
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J A<l>^  / 


(2) 


1/ 


as  a measure  of  fractional  time  or  frequency  stability  in  the  region  of 
measurement  times  from  about  1 to  10^  s where  the  effects  of  low 
frequency  divergence  can  often  be  ignored  [7,8].  The  squared  phase 
deviation  of  the  zero  crossings  is  given  approximately  by [8] 


s'"  (Sy  (f))/f"  df 

1/27t^  r 


(3) 


Figure  2.  Sy(f)  of  the 
probe  source  locked  a 
reference  resonance.  For  A 
G(f)  = (l/(10f) , for  B G(f) 
= (l/(10f))(l  + l/(40f)), 
for  C G(f)  = l/(10f))  (1  + 
l/(40f))(l  + l/(160f)). 
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Figure  3.  The  solid 
curves  show  (7y(r)  and  the 
dashed  curves  show  mod 

versus  r for  curves 
A,  B,  and  C of  Fig.  2. 

The  a,  b,  and  c points 
show  estimations  of  mod 
ay(r)  using  Eqs . 2 and  3 
with  fjj  = fho  • Also  shown 
is  CTy° (r)  and  Oy^ (r) . 


for  a sample  time  r and  noise  bandwidth  fj^  . For  a noise  spectrum  which 
varies  as  Sy(f)=  K'f"  where  n > 2,  the  integral  is  dominated  by  the  high 
frequency  bandwidth  even  for  very  long  measurement  times . The 

fractional- frequency  (or  time)  stability  given  by  Eq . 2 decreases  as 
just  as  does  Oy{r) , due  to  the  increase  in  measurement  time  and  not  to  a 
decrease  in  the  value  of  . The  integral  in  Eq.  3 can  be  integrated  by 
parts  for  the  various  segments  of  Sy(f).  This  makes  it  easier  to 
evaluate  and  optimize  the  performance  of  the  overall  system  than  by  using 
a process  based  on 

The  contribution  of  the  high  frequency  noise  can  be  greatly  reduced  by 
phase  averaging  the  data  points [9].  The  data  at  measurement  time  r = nr^ 
(where  is  the  data  interval)  is  obtained  by  averaging  the  n adjacent 
phase  points.  This  is  equivalent  to  using  mod  Oy(r)  to  analyze  the 
data [ 9 , 11 , 12  ] . The  standard  expression  for  mod  <7y(r)  contains  an 
enormous  number  of  terms  and  is  quite  laborious  to  compute [ 11- 12 ] . It 
[13-15]  has  been  pointed  out  that  mod  ay(r)  can  be  reduced  to 


mod  ay ( r ) 


2 


4 2 2 


Sy  (f)sin®  (Trr^nf) 

df. 

f^  sin^  (ttTq  f ) 


(4) 


which  is  much  more  manageable.  Nevertheless  it  still  requires  numerical 
calculations  to  determine  which  segment  of  the  phase  noise  dominates  the 
integral.  We  can  estimate  mod  ay(r)  from  Eqs.  2 and  3 by  using  fj^  = 
f^jo/n,  where  fj^^  is  the  hardware  bandwidth  of  the  measurement  system  and 
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T = nr^  is  the  measurement  time.  The  integral  for  in  Eq.  3 is  now 
substantially  reduced  for  measurement  times  large  enough  that  fj^^/n  is 
less  than  the  bandwidth  of  the  servo  system.  This  is  in  contrast  to  the 
original  calculation  for  Eq.  3 where  the  integral  must  always  increase 
with  T . This  integral  can  also  be  divided  into  parts  and  integrated 
analytically.  The  fractional  frequency  stability  computed  for  A,  B,  and 
C of  Fig.  2 using  mod  o^{t)  are  shown  as  the  dashed  curves  A' , B' , and  C' 
of  Fig.  3.  The  points  labeled  a,  b,  c show  the  results  of  estimating  mod 
ay(r)  using  Eqs . 2 and  3 with  f^^  = fj^^/n.  The  agreement  with  the  dashed 
curves  is  very  good.  The  calculations  for  and  mod  Oy{r)  are 
virtually  independent  of  using  an  upper  cutoff  frequency  for  the 
integration  or  a simple  low  pass  filter  of  the  same  bandwidth.  The  use 
of  mod  cTy(r)  reduces  by  a factor  of  about  100  the  time  necessary  to  reach 
the  performance  of  the  reference  plus  the  interrogation  noise,  which  in 
this  case  limits  the  performance  for  times  longer  than  100  s. 

3.  Discussion 


The  primary  utility  of  these  results  is  the  insight  into  the  origin  of 
the  major  contributions  to  Oy{r)  and  mod  £Jy(r)  as  a function  of  the  servo 
gain  and  the  analysis  bandwidth,  and  the  effect  of  narrowing  the 
bandwidth  with  longer  measurement  times.  The  specific  examples 
illustrate  that  it  is  generally  necessary  to  use  a second-order  loop  to 
lock  the  probe  source  to  the  reference  resonance,  but  that  a third-order 
loop  offers  little  additional  improvement  when  drift  in  the  probe  is  not 
serious.  In  addition,  if  the  reference  resonance  is  substantially  more 
stable  than  the  probe  source,  it  can  be  very  useful  to  use  mod  Oy{r)  to 
analyze  the  output  of  the  probe  source.  We  have  introduced  a simple 
measure  of  frequency  stability  that  is  easy  to  use  for  optimizing  the 
servo  and  the  analysis  system.  With  the  techniques  introduced  here  it  is 
possible  to  realize  fractional  frequency  stabilities  and  time  prediction 
of  the  local  oscillator  which  are  characteristic  of  the  reference 
resonance  in  a way  that  is  orders  of  magnitude  faster  than  those  using 
more  traditional  approaches. 
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